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Abstract. We describe how to apply the transport method to compute inflationary observ¬ 
ables in a broad range of multiple-field models. The method is efficient and encompasses 
scenarios with curved field-space metrics, violations of slow-roll conditions and turns of the 
trajectory in field space. It can be used for an arbitrary mass spectrum, including massive 
modes and models with quasi-single-field dynamics. 

In this note we focus on practical issues. It is accompanied by a Mathematica code which 
can be used to explore suitable models, or as a basis for further development. 
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1 Introduction 

Inflation [1—3] is a scenario for the very early universe according to which all large-scale struc¬ 
ture originated as quantum fluctuations. This idea has been tested by increasingly detailed 
observations of anisotropies in the cosmic microwave background (CMB), and its broad pre¬ 
dictions are now known to be compatible with their measured statistical properties [4]. This 
is a significant achievement which symbolises the maturation of modern cosmology into a pre¬ 
cision science. However, despite this phenomenological success, it remains unclear whether 
or not the microphysical origin of inflation can be understood. 

The simplest inflationary models comprise only a single scalar degree of freedom, taken 
to have a canonical kinetic term and a potential representing self-interactions. Minimal mod¬ 
els of this kind are sufficient to obtain predictions consistent with observational constraints, 
and therefore are natural from the perspective of simplicity and economy. From the perspec¬ 
tive of fundamental physics, however, their status is uncertain. If the inflationary scale is 
widely separated from the next relevant mass scale it may be natural to have only a single 
degree of freedom which is light by comparison, and a simple potential. But if the inflation¬ 
ary scale is not far below the next relevant scale—in particular, if the ultraviolet completion 
of the theory gives rise to multiple degrees of freedom which are light compared to the infla¬ 
tionary scale—then it may be that a model with several active scalar fields, coupled through 
a nontrivial potential or field-space metric, represents the most natural possibility. 

Which of these choices is a better match for the large-scale structure we observe in 
our universe is a question to be resolved by measurement. To do so will require a clear 
understanding of the qualitative predictions which can be obtained in each scenario. In 
single-field potential-dominated models, an extended campaign of exploration has provided 
guidance about what can be expected. In this case the textbook approach to perturbations 
is applicable. Only a few e-folds around horizon exit are relevant; the decoupling principle 
tells us that effects long before horizon exit are negligible, and long after it the perturbations 


- 1 - 




become conserved. During these few e-folds we can assume the Hubble rate to be nearly 
constant and take the inflaton mass to be negligible. 

In more complex models these simplifications no longer apply. If the field-space trajec¬ 
tory exhibits turns or other features around horizon exit then the calculation must usually 
begin far inside the horizon, tracking effects from perturbations with masses around the Hub¬ 
ble scale. Also, where multiple fields remain relevant after horizon exit we must continue to 
integrate until their observable effects are extinguished, possibly much later in the inflation¬ 
ary era or even long after it has ended. All but the heaviest modes will be relevant after 
horizon exit, including perturbations in the momenta. 

These complexities frustrate the traditional textbook approach. If only a few of them 
are present it may still be possible to make analytic progress. But more generally, in models 
where they all occur, a numerical method is almost essential. 

In this paper we show how the complexities introduced by many relevant scales can be 
accommodated, whether these scales are associated with masses in the particle spectrum or 
curvature scales in the field-space manifold. We illustrate our method with an explicit Mathe¬ 
matical implementation. 1 The version discussed here is available from transportmethod.com. 
We have attempted to simplify it as much as possible, with the intention of making it easy 
to follow. It can be used to compute the two-point function of inflationary fluctuations in a 
model with any number of fields and arbitrary potential and field-space metric. The numer¬ 
ical method (to be described in §2) is efficient, and therefore despite being implemented in 
Mathematica it is fast enough for practical model exploration. In situations where Mathemat- 
ica is too slow or inconvenient, such as Monte Carlo sampling, it could serve as a reference 
implementation. 2 

Synopsis. —This paper is divided into four principal parts. In §2 we explain how to derive 
differential equations which express the time evolution of field-space correlation functions 
in the spatially flat gauge. We allow an arbitrary potential and field-space metric. With 
appropriate initial conditions this system of equations can be used to capture contributions 
(including quantum effects) from all mass scales as the fluctuations approach, pass through, 
and eventually evolve outside the Hubble length. We discuss the selection of initial conditions 
in §3. 

In §4.1 we explain how to relate the flat-gauge field-space correlation functions to the 
statistical properties of the density perturbation, which is the observable quantity. 

For a given model the major numerical uncertainty is the duration of nontrivial evo¬ 
lution on super-Hubble scales. In principle—no matter which scheme we use to compute 
the properties of observables—the equations for all inflationary perturbations should be in¬ 
tegrated up to the last scattering surface, where they supply initial conditions for the cosmic 

1 Tliis implementation was originally developed to study a model of D-brane inflation [5]. The precise model 
is described in Ref. [6]. The Lagrangian is C = a 2 $ — V{(j>)\ where a is the scale factor and T 3 is the 

brane tension. It consists of six fields tj> % representing coordinates in the throat of a Klcbanov-Witten geometry 
which can be described by a noncompact conifold built over the five-dimensional SU(2) x SU(2))/U(1) coset 
space T 1 ’ 1 . The details of this geometry are encoded in the nontrivial field-space metric . The potential 
includes stochastic contributions from the bulk and it consists of ~ 600 terms. 

2 Other general purpose codes exist. Fieldlnf is a Fortran code capable of computing the inflationary power 
spectrum in an JV-field model with a nontrivial field-space metric [7-9]. ModeCode and MultiModeCode are 
similar Fortran codes designed (respectively) for single- and multiple-field models. They are restricted to 
a trivial field-space metric but emphasize Monte Carlo sampling [10-13]. Pyflation is a Python code which 
solves the Mukhanov-Sasaki equation to first-order in multiple-field models, and to second-order in single-field 
models [14-16]. Like ModeCode it is restricted to a trivial field-space metric. 
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microwave background anisotropies. In practice this is very onerous, and anyway would re¬ 
quire us to integrate through epochs of cosmological history, such as reheating, about which 
we know nothing. To evade both these issues we must usually rely on the dynamics becoming 
‘adiabatic’ at some point during inflation, or not long after—meaning that the isocurvature 
modes which can source time dependence of the density perturbation become exhausted, and 
it ceases to evolve. In §4.2 we discuss the issues which arise when trying to detect whether 
this limit has been reached. 

Notation. —We set c = h = 1 and express the gravitational coupling by the reduced Planck 
mass Mp 2 = 8irG. Greek indices from the beginning of the alphabet, (a, /3, ...) label the 
species of light scalar fields; Greek indices from the middle of the alphabet (//, u, ...) label 
spacetime dimensions. Spacetime indices are not needed except in Eq. (2.1). 

2 Transport equations for correlation functions 

An inflationary model with curved field space is governed by the action 

S= l -J d 3 xdt V=g{M$R - G a pg^d^ a dA 0 ~ 2V }, (2.1) 

where G a ^ is the field-space metric, V is the potential and g^ u is the space-time metric. At 
background level we take this to be Robertson-Walker, 

ds 2 = —dt 2 + a(t) 2 dx 2 , ( 2 . 2 ) 

where a(t ) is the scale factor. Both G a p and V may be arbitrary functions of the fields c/> a 
provided they are compatible with field configurations which realize an inflationary epoch. 
The equation of motion for the unperturbed background fields is 

V t 4> a + 3 H^ a + G afi Vp = 0, (2.3) 

where Vp = dpV, an overdot denotes partial differentiation with respect to cosmic time t and 
T>t denotes a covariant time derivative, VtX a = X a + T^^X'y. The connexion F^ is the 
Levi-Civita connexion compatible with G a p. 

Our aim is to study quantum fluctuations in this model, which are characterized by 
correlation functions of the independent degrees of freedom. Their precise identity is in¬ 
fluenced by our choice of spacetime coordinates. In this paper we define time t so that 
slices of constant t have zero Ricci curvature, up to possible tensor modes which we neglect. 
In these coordinates the independent degrees of freedom are fluctuations 5(f> a in the fields, 
and the correlation functions characterizing quantum fluctuations are the n-point functions 
(50 Q! (ki, ti)6(j)^(k 2 , £ 2 ))) (<5</>“(ki, t\)5(f)P(k- 2 , t 2 )^^ 7 (k 3 , £ 3 )), and so on, together with their 
derivatives. For applications to inflation we typically require only the equal-time case where 
all ti are evaluated at some common point t. The expectation value (• • •) is taken in a state 
which coincides with the Minkowski vacuum on deeply subhorizon scales. 

2.1 Capturing physical effects from all mass scales 

In order to ensure that we capture relevant physics from all mass scales, we begin the calcula¬ 
tion sufficiently far inside the horizon that vacuum initial conditions apply. As we will show 
below, in de Sitter space all degrees of freedom of fixed mass become effectively massless on 
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subhorizon scales, so—irrespective of the mass spectrum—we can obtain initial conditions 
for each n-point function to arbitrary accuracy by beginning the calculation sufficiently long 
before horizon exit. The details are discussed in §3. 

We then apply the in-in formalism to derive an evolution equation for each n-point func¬ 
tion, incorporating all masses and (in principle) quantum effects. 3 These evolution equations 
are equivalent to the separation into in-out expectation values and subsequent Feynman 
expansion used by Maldacena and later authors to obtain analytic estimates of the corre¬ 
lation functions [19]. But unlike the expansion into diagrams they do not involve Green’s 
functions—only ordinary differential equations. Therefore they constitute a differential for¬ 
mulation of the theory rather than an integral one, and can be handled by conventional ODE 
solvers. 


Perturbed action. —To second order in amplitude, the action governing small fluctuations 
5(j) a around a homogeneous solution (j) a (t ) of (2.1) can be written [20, 21] 


S ~IJ (03 dt « 3 {^[W“(k)] [V t 8<l>P{- k)] 

where the effective mass matrix M a a is defined by 


G a p + M a A5ct> a ( k)<l/(-k) , 


(2.4) 


M a p = V a .p - 


1 

a 3 Mp 


Vt 



(2.5) 


In this expression V Q -p = dpV a — T^Fy is the covariant derivative of V a , and R a \^p is the 
Riemann tensor built from the metric connexion TT . As before, an overdot denotes a partial 
derivative with respect to t. Both M a p and G a p should be evaluated on the homogeneous 
background 4> a (t). 

Our formalism requires only Eq. (2.4). It is not necessary that it derives from an 
action of the form (2.1) which controls both the background and fluctuations. It particular, 
it may happen that (2.4) applies to the fluctuations in scenarios which have a more general 
noncanonical kinetic structure than (2.1). (Note, however, that (2.4) is not sufficiently general 
to cover fluctuations in a P(X , cj)) model where the Lorentz symmetry between time- and 
space-derivative terms would be broken by a nontrivial sound speed Cg.) Where Eq. (2.4) 
applies, our evolution equations for the two-point correlation functions apply likewise. They 
may be used to compute the properties of the fluctuations, although the background equations 
[Eqs. (2.8)-(2.9) below] would require modification. 

The constituents of Eq. (2.4), including the perturbation 5(/) a , transform covariantly 
under a change of coordinates in field space. This implies that 5<j) a must be understood as 
a tangent vector, not a coordinate displacement. The necessary formalism underlying this 
interpretation was given by Gong Sz Tanaka [22]; see also Ref. [20]. To lowest order in 5(f> a 
this makes no difference, but it would become important in any attempt to extend (2.4) to 
third order or above. 


Quantization. —To quantize the fluctuations we define a momentum 5p a by the rule 5p a = 
5S/5(T> t 5(j) a ). Then 5(j) a and 5p a are to be treated as operators satisfying the commutation 

3 In principle the evolution equations contain the same information as the loop expansion of conventional 
perturbation theory, but not nonperturbative information such as instanton effects. In practice one must 
truncate each evolution equation, which is equivalent to truncating the loop expansion at a particular order. 
In this paper we work to tree level, which is already sufficient to capture those quantum interference effects 
around horizon crossing which determine the ‘quantum’ part of the Feynman calculation [17, 18]. 
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relation 


[ty a (ki),<S P/ j(k 2 )] = i(27r) 3 ^<5(k 1 + k 2 ) 


( 2 . 6 ) 


and the Hamiltonian is 

n = S W) 3 ([W a ( k )]^-(- k ) -c), (2.7) 

where C is the Lagrangian density appearing in (2.4). The equations of motion for 5cj) a and 
5p a can be determined from T-L via the Heisenberg equation. 

In practice it is numerically more convenient to integrate in terms of the e-folding 
number dlV = H dt rather than df itself. To do so at the level of the background we define 

d 6 a 

7T ° = w = VN(pa ’ (2 ' 8) 

in which it should be remembered that <p a (being a coordinate) behaves like a field-space 
scalar, whereas ir a (being the derivative of a coordinate) behaves like a field-space vector. 
The background equations of motion now comprise (2.8) together with an evolution equation 
for 7r“, 

G^Vr 

2W* = ( e -3)7r Q --^. (2.9) 

To effect a similar change for the quantized perturbations we define 

S7T a = ^ = V N S<j> a , (2.10) 

where the index on 5<j) a should be lowered using the metric G a p. Because the operations 
of taking covariant perturbations and covariant time-derivatives commute, it follows that 
5{T>N(j) a ) = T>n 5cj) a , and therefore 5n a can be regarded as an honest perturbation to the 
rescaled background field (2.8) [22], This identification is not spoiled by raising or lowering 
the index because the metric is covariantly constant. The equations of motion for 5(j) a and 
5n a can now be written 

V N 8<p a = - '[Sc^.U] (2.11a) 

11 

V N 6n a = ~^-[dTr a ,n] + (e - 3)<5 tt“. (2.11b) 

Similar operator equations will hold in the quantum theory, possibly modified by renormal¬ 
izations required to define composite operators appearing in the commutators [•,%]. Because 
these are operator equations they hold for any insertion of 5(j) a or 5n a in a correlation func¬ 
tion, provided it is not coincident with any other operator. If we work only to tree-level the 
complexities associated with renormalization are not needed, and we can work directly with 
the bare equations. The noncanonical term in the evolution equation for 5n a arises from the 
explicit time-dependent factors a 3 and H which appear in (2.10). 

Transport equations. —Salopek, Bond Sz Bardeen pointed out that a single solution of 
the 2 N differential equations (2.11) is not sufficient to compute the two-point correlation 
functions [23]. A single solution characterizes only how the late-time 5(j) a and 5ir a pertur¬ 
bations respond to a particular linear combination of fluctuations at an earlier time, and to 
compute a correlation function we must know how the late-tinre perturbations respond to 
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an arbitrary early-time perturbation. This entails calculating 2N solutions of (2.11), one for 
each independent initial condition. 

Various formalisms exist to compute the required solutions . 1 We choose to use the 
operator equations (2.11) to obtain evolution equations for each n-point function. There is 
no loss of information compared with solving for the held modes themselves, because only the 
correlation functions are meaningful: the predictions of an inflationary model are statistical, 
and are obtained by interpreting the correlation functions as ensemble averages and supposing 
that our particular universe is typical. 

In this paper we deal only with the equal-time two-point functions, which are sufficient 
to obtain lowest-order inflationary observables. There are four such functions: (<5</> a <5</U), 
(57r a 8(j)P), (8(j) a 5irP), and (d^STr 13 ). To compress notation we denote a generic perturbation 
such as 5(f> a or 5n a as X a . The index a ranges over the held and momentum perturbations 
for each species a. To distinguish these we continue to label held perturbations by a, /3, , 

but add a bar to the species label for a momentum perturbation, giving a, j3, ..., and so on. 
A generic equal-time two-point function can now be written 

y'ab 

(X a (k 1 )X\k 2 )) = (2vr) 3 5(ki + k,)-^-. (2.12) 

The evolution equation for T, ah follows from Ehrenfest’s theorem, 

V N (X a X b ) = (( V N X a )X b ) + ( X a (V N X b )). (2.13) 

This argument, valid for quantum-mechanical correlation functions, was given by Mulryne [32]. 
It entails solving of order 2 N x 2 N differential equations with a single initial condition, con¬ 
sistent with the counting argument given above. Using the symmetries of (X a X b ) shows 
that there are 3 N(N + l)/2 independent equations, assuming that we use the Weyl-ordered 
correlation functions to be defined in §3.1. 

Up to this point our treatment has been exact, except that by implicitly computing 
expectation values in a single state, corresponding to a hxed held configuration, we have 
restricted attention to what is visible in perturbation theory in the vicinity of that configu¬ 
ration. Therefore Eqs. (2.11) and (2.13) contain the same information as the loop expansion, 
and Eqs. (2.11) must contain information about all orders in interactions, represented by 
terms of all orders in 5(f> a and 5n a on the right-hand side. Hence 

V N X a = u a b X b + ■■■ , (2.14) 

where u a b is a matrix which can be computed 4 5 using Eqs. (2.11), and ’ denotes terms 
of order X a X b and higher which we have omitted. Neglecting these terms corresponds to 

4 Salopek, Bond & Bardeen decomposed the late-time fluctuations into a linear combination of creation- 
annihilation operators for the early-time fields, and solved for the resulting mixing matrix [23] . See Ringeval [9] , 
Huston & Christopherson [16] and Price, Frazer, Xu, Peiris & Easther [13] for recent applications. McAllister, 
Renaux-Petel and Xu solved Eqs. (2.11) explicitly, once for each independent initial condition [24]. Lalak, 
Langlois, Pokorski & Turzynski applied a method very similar to that proposed by Salopek et al. [25]. The 
T-matrix or ‘propagator’ introduced in Ref. [26] is of a similar kind. Rigopoulos, Shellard & van Tent [27, 28] 
elaborated a formalism due to Groot Nibbelink & van Tent [29], using a basis aligned with the instantaneous 
background trajectory. This can reduce the number of integrations required if we are prepared to give up 
knowledge of the isocurvature modes. A similar approach was used by Peterson & Tegmark [30, 31]. 

5 As explained above, we are neglecting renormalization and operator mixing effects which may be generated 
by the proper definition of composite operators beyond tree level. 
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working at tree-level in the loop expansion. On the left-hand side, the covariant derivative 
Vn acts on the generic label a appropriately for the ‘barred’ and ‘unbarred’ types, so that 
V N X a = d N X a + T^ir^X'r and V N X 5 = d N X & + T^tt^XT 
Combining Eqs. (2.12), (2.13) and (2.14) gives 

V N Z ab = u a J7 cb + u b c Y? c + • ■ ■ , (2.15) 


which we describe as the transport equation for T, ab [33, 34], The curved field-space version 
derived here was first given in Ref. [20]. Comparison with Eqs. (2.4), (2.7) and (2.11) shows 
that the components of satisfy 


u a p 


u 



0 

5 $ 

cq & 

P a 2 H 2 H 2 

(e - 3)8p. 


(2.16) 


Mass dependence. —The epoch of ‘horizon exit’ occurs when the physical wavelength of 
order a/k associated with the comoving wavenumber k becomes comparable to the Hubble 
length 1/H. At this time the ratio k/aH is of order unity. Prior to horizon exit k/aH > 1, 
and provided we are not too close to the start of the inflationary era it will be possible to 
find a point where (k/aH) 2 is much larger than any component of the effective mass matrix 
M a p/H 2 . If we choose to begin the calculation at or before this time then all fields can 
be treated as effectively massless. Where horizon exit is too close to the start of inflation 
it will not be possible to make M a p/H 2 entirely negligible, and we must find some other 
way to supply initial conditions, presumably depending on the pre-inflationary history. The 
calculation then becomes model dependent, but not harder as a matter of principle. In this 
paper we will not consider such possibilities. 

If H and M a p are nearly constant and the components of M a p are at most a few orders 
of magnitude larger than H 2 * * then the point where all fields become effectively massless might 
lie no more than N > 3 e-folds before horizon exit. At the other extreme, if H ~ 10 12 GeV 
(corresponding to roughly GUT-scale inflation) but M a p contains terms of order Mp then 
it could be necessary to begin N > 14 e-folds before horizon exit. These estimates require 
refinement if H or M a p vary significantly (see §3.2 for a numerical prescription). After hori¬ 
zon exit, k/aH becomes exponentially small and all but the most suppressed contributions 
to M a p will be relevant. 6 

Eqs. (2.15) and (2.16) provide a unified way to study both sub- and super-horizon 
regimes while retaining all relevant contributions to M a p. In the literature these regimes are 
sometimes associated with ‘quantum’ and ‘classical’ behaviour, but both of these descriptions 

6 In certain models there may be a superheavy scale H above which all modes can be neglected: the 

fluctuations in these modes decay rapidly because of their large mass. Also, the potential for a superheavy 

field is so steep that the background trajectory can be assumed to make no excursion in its direction. 

Although such a superheavy scale is normally assumed to exist, it has recently been appreciated that it is 

not straightforward to decide how large a mass is required before a field-space direction is negligible in this 

sense. The effect of massive modes is suppressed by inverse powers of the heavy mass M, but the rate of turn 

of the trajectory can be large. This gives a large number which can compensate for the smallness of 1/M, 
making the direction more relevant than it would appear. A literature has developed to study these effects; 
for example, see Refs. [35, 36]. 
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are marginally misleading. In the subhorizon era we work only to tree level and therefore 
true quantum effects are absent, but the initial conditions are quantum-mechanical and mix 
growing- and decaying-mode solutions for the elementary wavefunctions which contribute to 
Y, ab . It is the interference between these modes which determines the higher-order correlations 
imprinted around the time of horizon exit. In the superhorizon era the evolution becomes 
classical in the restricted sense that decaying solutions die away. 


‘In—in’ and i 5N i limits. —Once suitable initial conditions have been selected, it does not 
matter what spectrum of mass scales exists in M a p, or whether H varies significantly during 
or after the epoch of horizon exit. Eqs. (2.15)-(2.16) provide an alternative to the full 
diagrammatic description of the in-in formalism, but one which is equivalent. No further 
approximations are required. To determine the evolution of each correlation function we need 
only integrate the transport equation. 

When written in this form it is simple to obtain the connexion between the in-in for¬ 
malism and the ‘separate universe picture’, which gives an intuitive classical description of 
the evolving fluctuations on superhorizon scales [17, 21, 37-40]. In this limit the transport 
equation becomes a Jacobi equation describing the dispersion of neighbouring inflationary 
trajectories in held space and can be integrated analytically to produce the well-known l 5N’ 
Taylor expansion [20, 26, 32]. 7 


Scale dependence of two-point function. —A similar transport equation can be obtained 
for the scale dependence of the 2 -point correlation function, which we measure using the 
matrix 


n 


ab _ 


dE 


ab 


din A: 


(2.17) 


A transport equation for n ab can be obtained by differentiating Eq. (2.15) [42] 


V N n ab = 


din A; 


V N Yt ab = u° 


n cb + u b , 


n ac + 


du a c 
d In k 


Yf b + 


d u b C ] 

d In k 


(2.18) 


Tensor modes 8 .- Tensor perturbations 7 y are transverse, traveless perturbations of the 
spatial metric representing gravitational waves. Up to second order in amplitude their fluc¬ 
tuations are controlled by the action 

M 2 r f k 2 ) 

S D / d 3 xdt a 3 7 . 7 , - ^Hjlij | (2.19) 

where the Latin indices i, j run over the three spatial coordinates. To obtain scalar equations 
it is convenient to decompose 7 ^ into a basis of polarizations. In Fourier space this gives 

/ d 3 A- , 

( 2 . 20 ) 

where the polarization sum s runs over the orthogonal states s G { + , x}. The corresponding 
polarization matrices are traceless and satisfy A: i e| J (k) = 0, and are normalized so that 
e| jefj = 25 SS '. Each polarization Mp7 s (k)/\/2 behaves as a canonically-normalized free scalar 

'It is possible, but substantially more complex, to see how the l SN’ description emerges from the diagram¬ 
matic expansion [41]. 

s The calculations reported in this section were performed in collaboration with Sean Butchers. 
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field. Therefore the two-point correlation function of tensor perturbations is insensitive to 
the mass hierarchies of the system. 

To write an evolution equation for it we define the tensor momentum ir s = d 7 s /dlV and 
collect 7 s and i r s into a two-component vector 7 ® = ( 7 s , 7 r s ). The labels 0 , lb, ..., range over 
the tensor polarization and its momentum. We write the two-point function of 7 * as 

(7”(k 1 ) 7s 1 i(k 2 )) = (27r) 3 4 s ^(k + k')T*. (2.21) 


It follows directly from Eqs. (2.15) and (2.16) that r 0b obeys the transport equation 

d T ob 


d N 


= T 0C + • 


( 2 . 22 ) 


where (with no summation implied) = 0, = 1, w n j = — k 2 / ( a?H 2 ) and w^-n- = e — 3. 

Following the same procedure that lead to Eq. (2.18), it is straightforward to compute 
the scale dependence of the tensor spectrum. Defining the quantity 


„ <nb — 
71t = 


dT ob 
d In k ’ 


it follows that its equation of motion is 

dn T 0b 


_ „.,<u „ cb 1 ...id ojo; | 

— w c nx + w c nx + 


dN 


0c 1 pcb 


din k 


+ 


dw b , 


c p0C | 


din A; 


(2.23) 


(2.24) 
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2.2 Mathematica implementation 

In these grey panels we discuss the numerical Mathematica implementation of the transport 
method, available from transportmethod.com. In the control panel of this notebook one can 
specify the model to be evaluated and select which computations to perform. 

Which observables? —By default, the notebook computes observables at a chosen scale 
k. Computing the spectral index at this scale using Eq. (2.18) on subhorizon scales requires 
delicate cancelations between oscillating terms. For some models it can be considerably slower 
than computing E a& alone. In these cases it can be preferable to evaluate £ ob a few times (two 
is typically sufficient) and compute the spectral index via finite difference . a 

Power spectrum. —From the control panel one can also specify a range of fc-scales at which 
to compute £ ab , in order to obtain the power spectrum as a function of scale P^(fc). 

Example model. —Throughout this paper we use the 3-field model ‘number 2’ as an example. 
This is a simple extension of the case studied in Ref. [43] , which is a model of quasi-single-field 
inflation giving rise to a feature in P^(k) as a result of excitation of a heavy field via a nontrivial 
metric. (See §3 for a brief explanation of quasi-single-field scenarios.) 

Number 2 can be regarded as a ‘quasi-two-field’ example. We add another light field to obtain, 
in addition to the nontrivial behaviour coming from G a P , superhorizon evolution via a turn 
in the plane of the two light directions. This turn arrises in the region of field space where 
the metric is approximately the unit matrix. The turn occurs due to the hierarchy in the 
masses of the displaced fields — as the heavier field approaches its minimum, the direction 
of steepest descent becomes progressively more aligned with the lightest field. There is no 
direct motivation for this model, but at a qualitative level similar characteristics can arise in 
supergravity. Here we merely employ this example for illustrative purposes. 

The model has an equation of motion of the form (2.3). The potential is 

1 3 

v = - 2 Y.<<& ( 2 - 25 ) 

Z a-1 


and the mass ratios are m\/m\ = 30 and m 2 /m\ = 1/81. The metric takes the form 

(i r o\ 

G afi = T 1 0 , (2.26) 

VO 0 i ) 


with 



(2.27) 


“For practical purposes, however, it can be useful to solve Eq. (2.18) in order to understand how 
far inside the horizon we should start our computation; this will be discussed in §3.2. It is possible to 
select which method is used from the control panel. 


3 Initial conditions 

In a simple model we would typically apply (2.15) only in the superhorizon regime, where all 
masses are relevant because k/aH is exponentially small. We would estimate an initial value 
of £ ab for modes which are ‘light’ in the sense that k/aH dominates M a p/H 2 around horizon 
exit, and set all correlation functions to zero for ‘heavy’ modes to which this does not apply. 
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The justification is two-fold. First, heavy modes are orthogonal to the inflationary trajectory, 
so if this remains nearly straight throughout the epoch of horizon exit then fluctuations in 
these directions have no physical effect. Second, quantum fluctuations in massive modes 
decay exponentially, so they become irrelevant almost immediately. 

Under certain circumstances this approximation may miss effects from modes with in¬ 
termediate masses of order the Hubble scale or slightly larger. If bending of the inflationary 
trajectory is not negligible during horizon exit then fluctuations in massive modes can be 
partially converted into the adiabatic density perturbation before they have time to decay. 
Chen & Wang called this scenario ‘quasi-single field inflation’ [43-53]. It yields a distinctive 
bispectrum. Even if intermediate-mass modes are not relevant, turns in the field-space trajec¬ 
tory will cross-correlate fluctuations in different species. This can have a significant impact 
on the later evolution of observables. A striking example is the ‘destructive interference’ 
observed by McAllister, Renaux-Petel & Xu in Ref. [24]. 

As described in §2.1, we should account for all these effects by setting initial conditions 
sufficiently early that M a p/H 2 is negligible compared to ( k/aH ) 2 . Because the transport 
equation requires no approximation to be made regarding M a p the subsequent evolution is 
exact (at tree level) and will capture all quasi-single field and cross-correlation effects. 

In this section we explain how initial conditions can be computed in the deeply subhori¬ 
zon regime, where all correlation functions are dominated by kinetic contributions. 

3.1 Deep inside the horizon 

In this section we work in conformal time, defined by r = dt'/a(t'). 

Field correlation function. —In the massless limit, the Feynman two-point function eval¬ 
uated in the vacuum state is 9 [20] 

(S</> a (k',T')S^(k,r)) = (27r) 3 S(k+k') na H(T)H(r')(l-ikT)(l+ikr')e ik ( T - T, \ (3.1) 

where we have assumed r < r * & 7 . The quantity n q/ 3 (t',t) is the parallel propagator evaluated 
on the inflationary trajectory, 

n-V.r) = Pexp (- jT dr" T^{t")] ^ G^(t). (3.2) 

It transforms as a bitensor; the index a transforms as a tensor in the tangent space at (^(t 7 ), 
whereas the index /3 transforms as a tensor in the tangent space at The symbol 

P denotes path-ordering and rewrites its argument in order of decreasing time along the 
trajectory. 

We take the equal-time limit r = r 7 in which the parallel propagator reduces to the 
metric, and evaluate the correlation function well inside the horizon where \ k/aH\ ~ |fcr| 1. 
(If r represents a time N e-folds prior to horizon exit for the mode k. then \kr\ ~ e N .) This 
yields 

Tr2(^a/3 

{8ct) ol (k)8(j/{k!)) T » (27r) 3 h(k + k 7 ) ^ \kr\ 2 (3.3) 

9 The equal-time two-point function in a model with nontrivial field-space metric was calculated by Sasaki 

& Stewart [21]. The factor of the parallel propagator appearing in the unequal-time propagator was given in 

Ref. [20]. 
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in which H and G a @ should be evaluated at the common time r. 


Vacuum state. —Eq. (3.3) gives a suitable initial condition for all field-held correlation 
functions at sufficiently early times, but only if we assume that the mode of wavenumber k 
was practically in its vacuum state for at least a few e-folds before the time r so that vacuum 
initial conditions were applicable. 

This is not guaranteed. If inflation is sufficiently prolonged, a mode with fixed comoving 
wavenumber k must have originated on very small physical scales. Physics at these scales is 
presumably not governed by the effective theory used to describe inflation, so these short- 
scale modes will join it only when they are redshifted within its purview. Their state at that 
time should properly be regarded as a boundary condition needed to define the effective held 
theory. Like all details of ultraviolet physics, it cannot be predicted from within the effective 
theory. 

The inhuence of such boundary conditions was studied by Anderson, Molina-Paris & 
Mottola [54]. They found that, to be consistent with Einstein gravity as a low-energy de¬ 
scription, the effective stress tensor generated by unobserved high-energy fluctuations should 
correspond to sufficiently depopulated occupation numbers at large frequencies. If these oc¬ 
cupation numbers are conserved during redshifting then modes with these frequencies would 
join the effective description while practically in their vacuum state. In that case Eq. (3.3) 
will apply. This is the default assumption in many inhationary models. 

An alternative, studied by a number of authors (see eg. Refs. [54, 55]), is that modes join 
the effective description at a fixed time before horizon exit with nonzero occupation number. 
In this case Eq. (3.3) would require corrections. Whether this occurs is a model-dependent 
question, but if corrections are necessary they can be accommodated as a change of initial 
conditions. The transport equations themselves do not require modification. 

Correlation functions with momenta. —To compute correlation functions involving mo¬ 
menta we use the relation dIV = —dr/r, which implies TV = —tT> t . Differentiating the 
unequal-time field-held correlation function (3.1) and using that the parallel propagator is 
covariantly constant, 

TVII a/ V, T ) = £>rl l aP {r', t) = 0, (3.4) 


we obtain the unequal-time field-momentum correlation functions, 


(^“(k 7 , r')5(j)P( k, t)) 
(5cj) a ( k', r / )(57r /3 (k, r)) 


— (2-7r) 3 <5(k + H{r)H{r')\kT'\ 2 {l - i kT)e ik( - T ~ T '^ (3.5) 

2k 6 

— (2-7r) 3 <5(k + k 0 T) H{T)H{r')\kT\ 2 {\ + ikT)e ik{ - T ~ T '\ (3.6) 


We have neglected terms suppressed by the slow-roll parameter e = —H/H 2 , which are gen¬ 
erated by differentiation of H. In (3.1) we did not retain slow-roll suppressed contributions 
(although this can be done), so retaining them here would give an inconsistent set of cor¬ 
rections. Also, we have again assumed r < t' . The case r > t' can be obtained from these 
expressions by complex-conjugation. 

The equal-time limit can be obtained as before, but unlike the field-held correlation 
function the result is complex. However, the imaginary part vanishes at late times and 
therefore does not affect observables. For simplicity we can work with the symmetrized 
(or ‘Weyl ordered’) correlation function which is always real and coincides with the other 
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field-momentum two-point functions outside the horizon, 


^^(k)5/(k0 + j/(k')<57r Q (k)> T = -(27r) 3 «5(k + k0^^|fcr| 2 . (3.7) 

By a very similar procedure, the momentum-momentum correlation function can be found 
(also to leading order in slow-roll terms) to be 


Tj2/^iaf3 

(<57r"(k)<57^ /3 (k , )) = (2-7r) 3 <5(k + k 7 )——|fcr| 4 . 


(3.8) 


Collecting these results yields the universal initial conditions 


_ rr2/~fa/3 

= * * |fcr,[ 2 , Ef = Ef = - Zf = 


HlGf 


\krf, 


n. 


2 r ~'~' 7 “* 2 r "~' 7 “* 2 r “'~' 7 (3.9) 

f = HlGf \krf, nf = nf = -HlGf\krf, nf = 2HlG a f\krf 


where a subscript V denotes evaluation at the initial time. 

As tensor modes behave like free scalar fields (apart from a change in normalization) 
their initial conditions follow immediately, 


TT 2 

p77 — _ 51 1 U T 1 2 -p7T7 _ p7 7r - 

1 * j 2 I ™ 1 * I > 1 * 1 * 

77 __ 2 Hi 2 

n T* ~ nrl \^ T *\ > 








7 T'Y 777 

n T* = n T* = 


2 Hi 
' Ml 


| fcr* 


™T* = 


4 Hi 

Ml 


(3.10) 


| kr* 


Slow-roll corrections.—In principle, the initial conditions (3.9) should be corrected by 
slow-roll terms proportional to powers of e or its derivatives. Therefore although the transport 
equation (2.15) makes no use of the slow-roll approximation, our use of (3.9) does require 
that slow-roll is a fair approximation near the initial time. We do not need any form of 
the slow-roll approximation thereafter; the slow-roll conditions may be badly violated or fail 
entirely. 

In certain cases it may happen that a solution with initial conditions chosen to satisfy 
Eq. (3.9) will still relax to the correct solution, even if slow-roll is only marginally valid or 
weakly violated near the initial time, provided we begin the calculation sufficiently far before 
horizon exit. As for any numerical solution, some care may be required to check that results 
are stable to changes in the grid and the initial time. In practical calculations this implies 
that the initial time should be chosen so that it is comfortably earlier than any interesting 
dynamical effects which we hope to capture. 
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3.2 Mathematica implementation 

Increasing the number of e-folds of subhorizon evolution slows down the solver; but not starting 
the calculation sufficiently early leads to inaccurate results. The amount of subhorizon evolu¬ 
tion (as well as the accuracy settings of the Mathematica function NDSolve) required to obtain 
a sufficiently accurate result is model dependent, however it is useful to have a prescription 
that works as a guideline. 

One way to approach this is to use the behaviour of the spectral index. Sufficiently far inside the 
horizon, the dominant contribution to the scale dependence comes from Eq ' 3 oc k 4 . Therefore 
we expect the spectral index n s — 1 to be roughly 4 at early times, independent of the model. 
If the system is not given enough e-folds of subhorizon evolution, the delicate cancelations that 
occur for n a b will quickly break and give n s — 1 ^ 4. We have found that ensuring n s — 1 = 4 
for a sustained period of order ~ 3 e-folds is usually sufficient to obtain consistently good 
results. This method should not be used as a replacement for testing for convergence against 
changes in initial time and grid, but it can be used as an indicator. 


4 Connection with observables 


4.1 From field space to Q 

To calculate observables, the flat-gauge field and momentum correlation functions must be 
converted to correlation functions of the uniform-density gauge curvature perturbation (. 
Assuming all isocurvature modes decay, it is the curvature perturbation which sets initial 
conditions for density fluctuations in the later universe. In this section we briefly explain 
how an appropriate transformation can be extracted from the separate universe assumption 
using the methods of Ref. [56]. As explained there, the gauge transformation could equally 
well be derived from traditional perturbation theory in the large-scale limit. 

Gauge transformation. —According to the construction of Ref. [56], the number of e-folds 
AN between a point p on a spatially flat hypersurface at which the density is p p and an 
equivalent point on a nearby uniform density hypersurface with density p* is 


AN = 


d N 
dp 


( P* ~ Pp) + 


(4.1) 


where the omitted terms are higher order in p* — p p . Under a variation of p, it follows that 


i(AiV) » - A 


. dN 

, dp 


dp 


d(f> c 




dp 


dir 0 


Tp) + ■ ■ ■ 


(4.2) 


= N a 5ct) a + Nadir a + 


where ■ ■ ’ denotes terms of higher order in 5(j) a or Sir a which are not needed for the first- 
order gauge transformation. The last equality should be interpreted as a definition of N a 
and Na- No use is being made of the slow-roll approximation so the density p contains both 
potential and kinetic contributions. Performing the partial derivatives, we find 


N a = 
Na = 


2e V 

1 TTq 
2e(3 - e) Mp ' 


(4.3a) 

(4.3b) 


- 14 - 












The variation 5(AN) gives the fluctuation in e-folds required to reach p*, and therefore must 
be the curvature perturbation £. It follows that Eq. (4.2) expresses the gauge transformation 
to C at linear order. 

This argument was given in Ref. [20] to lowest order in the slow-roll approximation 
where the contribution for 5ir a can be neglected. In the formulation given here, Eqs. (4.2) 
and (4.3a)-(4.3a) apply to all orders in the slow-roll expansion. The argument of Ref. [56], 
generalized to a nontrivial field-space metric, shows that if desired we can use the Hamiltonian 
constraint to eliminate the 5iT a term. This would give 


2e Ml 


Na = 0 . 


(4.4a) 

(4.4b) 


Like (4.3a)—(4.3b), Eqs. (4.4a)-(4.4b) do not invoke the slow-roll approximation. 

Power spectrum. —The quantity of principal interest is the power spectrum, which is 
defined in terms of the equal time two-point function, 


(C(ki)C(k 2 )) = (27r) 3 d(ki + k 2 )^|. 


(4.5) 


In terms of the flat-gauge correlation functions, it follows from Eq. (4.2) that P ,£ can be 
written 


P c = N a N b T! 


ab 


(4.6) 


The L-dependence of the power spectrum can be handled similarly. We define the scalar 
spectral index to satisfy 

d In P c 


n, - 1 = 


din A: 


(4.7) 


k=k* 


where L* is the pivot scale. (The subscript St’ representing evaluation at the pivot scale should 
not be confused with the subscript V denoting evaluation at the initial time in Eqs. (3.9).) 
We find 

, N a N b d?, ab N a N b n ab 

- 1 = -fTdhU = Nj^- <4 ' 8) 

The running of the spectral index is defined to be 

dn s 


a = 


d In k 


(4.9) 


k=k * 


It can be computed using a finite-difference approximation, although another transport equa¬ 
tion could be written for it if desired. 


Tensor fraction. —The tensor power spectrum is defined by analogy with the scalar power 
spectrum 

(7ij(ki)7t?(k 2 )) = (2vr) 3 h(ki + k 2 )||. (4.10) 

The power in each polarization adds incoherently. Using the normalization condition e|-e®j = 
2 5 SS ' it follows that the total tensor power satisfies 

(7ij(ki)7»i(k 2 )) = X^S^7s(ki)7 s (k 2 ))e-jef_J = 2^(7 s (ki) 7s (k 2 )). (4.11) 

s s’ s 
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Each polarization is the same and therefore the final result is 4(7_|_7_|_), or equivalently 
4(7 x7x)- The tensor-to-scalar ratio r is defined to be 


and the tensor spectral index is 


r = 


P, 

Pc 


4P77 


kt 


d In Eh 


d In k 


k=k* 



p 77 * 


(4.12) 


(4.13) 


4.2 End of inflation or beyond? 

To complete the calculation, we must decide when to terminate the integration and measure 
final values for each observable. 

In principle we should track the evolution of all correlation functions up to the last 
scattering surface—or until just before horizon re-entry if we wish to study the assembly 
of large-scale structure. In practice this is extremely challenging, partly because the results 
depend on the unknown details of reheating and partly because we have little direct knowledge 
of the epochs between reheating and horizon re-entry. (For recent literature studying the 
impact of reheating on inflationary observables, see Refs. [57-60]. Ref. [61] is a recent review 
of reheating in general.) 

It is possible for £ to evolve whenever power remains in any ‘isocurvature’ modes, by 
which we mean phase space directions transverse to the inflationary trajectory. Microwave 
background data now strongly constrain the presence of isocurvature modes around the time 
of photon decoupling at z ~ 1100, but this provides only a lower limit on the decay time. We 
would normally aim to terminate the integration as early as is safe, to avoid being obliged to 
integrate through periods of the universe’s history where we must make assumptions about 
its evolution. This means we must be able to determine when all power in isocurvature modes 
has become exhausted—what is called the ‘adiabatic limit’ [62, 63]. 


Isocurvature power during inflation. —First consider evolution during the inflationary 
era. One option would be to track the two-point functions of each relevant degree of freedom, 
using Gram-Schmidt orthogonalization to construct linear combinations which measure the 
power transverse to the inflationary trajectory. For practical purposes we could suppose that 
the system is close enough to an adiabatic limit whenever all of these two-point functions 
become sufficiently small. 

This approach is feasible, but rather cumbersome. At least during slow-roll evolution 
an alternative is to use the optical analogy of trajectories flowing over field-space developed 
in Ref. [26]. When slow-roll is a good approximation there is no need to track momen¬ 
tum perturbations, and we can write an equation analogous to (2.14) purely for the field 
fluctuations, 

V N 5(j) a = w a p6(j)P, (4.14) 

where w a p is likewise an analogue of the expansion tensor (2.16), 

1 , / V a \ V 7 V X 

w a p = Vp(V N r) - X pV N ^V N ck x = Vp J - ( 4 -15) 

Here, V a = G a ^Vp. As explained in Ref. [26], the eigenvalues of w a p can be used to detect 
the presence of growing and decaying modes: a positive eigenvalue indicates a growing mode, 
whereas a negative eigenvalue indicates a decaying one. 
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Figure 1. Evolution of observables as a function of e-folds N: clockwise from top left, the power 
spectrum P^; the spectral index n s ; the running of spectral index a; and the tensor-to-scalar ratio r. 


In most circumstances one mode is constant or slowly evolving, and therefore gives an 
eigenvalue which is zero or slightly positive. Therefore, in an IV-field system, approach to the 
adiabatic limit is signalled by the appearance of N—l large negative eigenvalues. This method 
is simpler than computing all combinations of isocurvature correlation functions, but clearly 
shares its arbitrariness in deciding when the fluctuations have decayed sufficiently to declare 
that an adiabatic limit has been reached. There is always the possibility that extremely 
violent future dynamics could amplify even very small modes. 

After inflation. —This test applies only during slow-roll evolution, and therefore will typ¬ 
ically become unreliable some time before the end of inflation unless this is mediated by a 
sudden event such as a waterfall transition. When it applies, however, it may provide a 
rationale for terminating the integration at or before the end of the slow-roll phase. This is 
the best possible outcome. 

Much less can be said if slow-roll breaks down before complete decay of the isocurva¬ 
ture modes. In this case, one should follow the decay of the scalar species relevant during 
inflation into reheating products. Isocurvature modes may be transferred or amplified dur¬ 
ing this process. One must then begin a second integration, following the evolution of these 
fluctuations using suitable phase space coordinates; normally, the scalar species supporting 
inflation will no longer be the relevant variables, and the transport equation for their cor¬ 
relation functions will need to be replaced. The range of phenomenology which can occur 
during this post-inflationary phase is comparatively under-explored, and almost certainly 
model-dependent. 
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Figure 2. The power spectrum Pq as a function of e-folds N and wavenumber k. The range of scales 
corresponds to about 5 e-folds. Note the step in Pq(N) corresponding to superhorizon evolution, and 
the oscillatory features in P^(k) as a result of excitations of the heavy mode around horizon exit. 


4.3 Mathematica implementation 

As above, we use the model Number 2 as an illustration. 

Time dependence. —In Fig. 1 we plot the time evolution of the power spectrum, spectral 
index, running of the spectral index and tensor-to-scalar ratio for the scale leaving the horizon 
55 e-folds before the end of inflation. 

A turn in field space occurs after TV ~ 30 e-folds. The turn causes isocurvature modes to 
source evolution of P^, visible here as a large step. This corresponds to a sudden jump in the 
spectral index and running, and a corresponding drop in the value of the tensor-to-scalar ratio. 
Once the trajectory settles into the valley of the potential, no further evolution occurs. 

Scale dependence. —In Fig. 2 we plot the time and wavenumber evolution of the function 
P^(fc). Evidently each scale undergoes qualitatively similar evolution to the example shown 
in Fig. 1. There is a significant step around N ~ 30. However, it is also possible to see the 
effect of excitation of the heavy mode: it induces oscillatory features in P^{k) as a function of 
wavenumber. 

Approach to the adiabatic limit. —We plot the time evolution of the eigenvalues of the 
slow-roll expansion tensor w a p in Fig. 3. Around horizon crossing the metric is designed to 
excite the heaviest field, leading to a sudden enhancement of the isocurvature modes, here 
visible as a sharp excursion to positive eigenvalues. 

Soon after, one of the eigenvalues (the yellow line) becomes much more negative than the other 
two, indicating that isocurvature fluctuations associated with the heaviest field are decaying 
exponentially. The suppression is so rapid that the system is subsequently well-approximated 
by a two-field model. 

After 20 e-folds of superhorizon evolution the trajectory turns, causing excitation of the re¬ 
maining isocurvature mode. Its power is subsequently transferred to the adiabatic direction 
(approximately represented by the blue line), generating the step-like features in Fig. 1. After 
the turn the remaining isocurvature mode (the green line) is rapidly suppressed. At this point 
an adiabatic limit has been reached: the system has become effectively single-field, and £ is 
conserved. .. 0 








Figure 3. Eigenvalues of the slow-roll u a p as a function of number of e-folds. The left plot shows 
the whole evolution from horizon exit to the end of inflation. The right plot is zoomed in order to 
understand the transfer of power around the turn in field space. 

5 Final summary 

Computing precise predictions for observables in complex inflationary models requires tools 
going beyond textbook methods. This is especially true for models with a nontrivial field- 
space metric. 

First, it is helpful (but not mandatory) to use a covariant description of the system, 
especially when computing correlation functions of higher-order (see, eg., Ref. [20]). Second, 
to track the influence of curvature scales associated with the metric we require a computa¬ 
tional method which retains information about all mass scales in the problem. In this paper 
we describe a simple method for doing so, beginning with universal massless initial condi¬ 
tions long before horizon exit and solving a transport equation for the subsequent evolution. 
We have focused on the equal-time two-point functions because only these are required for 
the simplest inflationary observables. Nevertheless, only straightforward modifications are 
needed to compute unequal-time correlation functions or those of higher order. 

This method is applicable to a large class of models, including models descending from 
ideas in string theory or supergravity where nontrivial field-space metrics are often associated 
with a nontrivial Kahler potential. It gives a simple description which allows the analysis to 
proceed from deeply subhorizon scales to the end of inflation (or beyond), with no matching 
required around the time of horizon exit. The evolutionary equation does not make use of the 
slow-roll approximation—although our analytic initial conditions do—and therefore accounts 
for effects from violation of the slow-roll conditions, turns in field-space at any point on the 
inflationary trajectory, and the influence of massive fields or quasi-single-field dynamics. It 
does not yet apply to models with entirely arbitrary kinetic terms, such as P(X) or P(X, 4>) 
models. We leave these for future work. 

In this paper we have tried to explain how the transport method can be applied in 
practice to obtain predictions from inflationary models which exhibit one or more of these 
complexities. In addition we have attempted to highlight those points where our implemen¬ 
tation goes beyond standard textbook methods such as the separate-universe approximation, 
and also those situations to which the method cannot yet be applied. In these cases we 
have indicated whether the obstruction is a matter of principle, or just an artefact of current 
technology. Since our focus is on practical usage, we have provided a complete Mathematica 
implementation which is used as an example. The code is intended to be accessible. We hope 
it will serve both as a useful tool for investigating realistic models, and a platform to extend 
the range of scenarios for which predictions can be obtained. 
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